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Abstract. We apply general variational techniques to the problem of the counterion distribution around 
highly charged objects where strong condensation of counterions takes place. Within a field-theoretic formu- 
lation using a fluctuating electrostatic potential, the concept of surface-charge renormalization is recovered 
within a simple one-parameter variational procedure. As a test, we reproduce the Poisson-Boltzmann sur- 
face potential for a single charge planar surface both in the weak-charge and strong-charge regime. We 
then apply our techniques to non-planar geometries where closed-form solutions of the non-linear Poisson- 
Boltzmann equation are not available. In the cylindrical case, the Manning charge renormalization result is 
obtained in the limit of vanishing salt concentration. However, for intermediate salt concentrations a slow 
crossover to the non-charge-renormalized regime (at high salt) is found with a quasi-power-law behavior 
which helps to understand conflicting experimental and theoretical results for the electrostatic persistence 
length of polyelectrolytes. In the spherical geometry charge renormalization is only found at intermediate 
salt concentrations. 

PACS. 82.70.-y Disperse systems; complex fluids - 61.20.Ja Computer simulation of liquid structure - 
61.20.Qg Structure of associated liquids: electrolytes, molten salts, etc. 



1 Introduction 



The behavior of charged systems has attracted renewed 
interest in the last few years H. One of the topics of in- 
terest is centered around the effects of multi-valent coun- 
terions, in which case correlations between ions become 
important and mean-field type theories break down. In 
such situations it is known experimentally and theoreti- 
cally that equally charged bodies can attract each other. 
But even for monovalent ions and moderately to highly 
charged surfaces many open questions remain: In such sit- 
uations the coupling parameter (which measures to which 
extent corrections to the mean-field or saddle-point so- 
lution are important) is not very high, but the surface 
potential can exceed thermal energy by far such that non- 
linear effects become important. The complication in this 
case is that the mean-field or Poisson-Boltzmann equation 
is a non-linear differential equation which in the presence 
of salt ions can be solved in closed form only in the planar 
geometry In the cylindrical and spherical geometry a 
number of different approximations have been proposed, 
which all more or less agree on the fact that counteri- 
ons condense in the vicinity of the charged surface such 
as to reduce its effective surface charge Q^]. The original 
Manning condensation argument states that the effective 
charge of a charged cylinder is maintained at a constant 
level of one charge per Bjerrum length^]. The prescrip- 
tions for calculating the effective charge are numerous and 



so are the predictions for the precise value of the effective 
charge as a function of external parameters such as salt, 
temperature, pH and so on[[3|[7],||]. 

In this paper we show how to approach the problem 
using a field-theoretic formulation of a variational theory, 
within which the effective surface charge is treated as a 
variational parameter. There are in principle two different 
approaches, one is based on the standard density func- 
tional theory, within which the counterion distribution 
is variationally determined. This approach is clumsy for 
more complicated geometries and in principle is difficult 
to generalize such as to include strong-coupling (i.e. de- 
viations from mean-field) effects. We therefore developed 
an alternative approach based on an exact field-theoretic 
formulation using the fluctuating imaginary electric po- 
tential. The surface charge is used as a variational pa- 
rameter and the variational field-theoretic action is Gaus- 
sian. This is probably the most unambiguous definition 
of the renormalized or effective surface charge, since no 
additional assumptions (other than that the variational 
Hamiltonian is Gaussian) is used. As a trivial byprod- 
uct, the effective interaction between charges is at large 
distance of the Debye-Hiickel type with renormalized ef- 
fective charges. Our scheme in principle allows to go be- 
yond mean-field, but we do not pursue such studies in this 
paper. As we show in the appendix, the standard Gibbs 
variational principle can be reformulated as a perturbative 
variational scheme which can be systematically improved 
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by a perturbation series (which we also leave as a prob- 
lem for the future). We test our approach for the planar 
charged wall, for which we essentially exactly reproduce 
the Poisson-Boltzmann result. We then tackle the charged 
cylinder. The effective cylinder charge interpolates contin- 
uously between the unrenormalized limit (at high salt con- 
centrations) and the Manning charge-renormalized limit 
(at low salt concentrations), in agreement with previous 
results 0]. However, this crossover occurs for fully charged 
polyelectrolytes over 4 orders of magnitude of the salt con- 
centration, and therefore it is important to take this slow 
crossover into account. For the controversial debate on 
the electrostatic persistence length it means that the stan- 
dard Odijk-Skolnick-Fixman result [p|[lO|,pd||, obtained on 
the linear (Dcbye-Huckel) level, is modified by the salt de- 
pendent effective polymer charge, which changes the salt 
dependency from Iqsf ~ K ~ 2 to £<jsf ~ k~ 1a , closer to 
experimental values]!^]. For the charged sphere the re- 
sultant behavior is different: Here charge renormalization 
takes place only at intermediate salt concentrations and 
both for high and low salt concentrations the standard 
Debye-Hiickel potential with unrenormalized charges is 
valid [pj. Our method can be easily generalized to more 
complicated geometries and situations. 



2 The model 

The partition function of a system consisting of N + posi- 
tively charged and iV_ negatively charged ions of the same 
valency q, interacting with an arbitrary charge distribu- 
tion <r(r), can be written as 



1 N+ 

7v+! n 

+ 3=1 



dr+ 



3 m) 



i n - 



3 = 1 



A? 



exp j- Y / drdr 'b/ 3 +( r ) _ <?/M r ) + cr(r)]«(r - r') 



[<lP+( r ') - IP- 



air')] 



q 2 £ B (N + + N_) 



v(0) (1 



The pair potential v(r) denotes the Coulomb interaction, 
v(r) = 1/r, but the first part of our discussion is rather 
general and independent of the specific form of v(r). The 
self-energy term in Eq.(|l|), proportional to g 2 ^_ev(0)/2, 
simply subtracts the diagonal term in the double integral. 
It will be later shown to renormalize the chemical po- 
tential (or fugacity) in a rather trivial way. The Bjerrum 
length Ib = e 2 / A-neksT is the distance at which two unit 
charges interact with thermal energy ksT. The length X t 
is the thermal wavelength, which results from integrating 
out kinetic degrees of freedom, the precise value of which 
is unimportant. The function S7(r) takes into account the 
presence of hard walls and restricts the positions of ions 
to an appropriate region in space, as will become more 
clear when we consider explicit examples in the following 
sections. The counterion density operators are 

N+ 



for positively charged ions with a similar expression for the 
negative ions. At this point, there is not much one can do 
with the partition function, mainly because the density 
operators p + and p_ enter in a quadratic fashion. The 
main step in deriving a field theory consists in introducing 
a unit operator which couples to the density operator p + , 

1 = fv P+ 5( P+ - p+) 



Vp+V^+ exp it / dr^+(r)[p+(r) - p + (r)} \ (2) 



(and a similar expression for the density of negative ions), 
where we used the integral representation of the delta 
function. This unit operator can be introduced into the 
partition function. Once introduced, it allows us to replace 
density operators /5 + (r) and p- (r) by the corresponding 
fluctuating density fields p + (r) and p_ (r) and the configu- 
rational integrals over the ion positions can be performed. 
As a result, the partition function now reads 



Vp + Vp_Vip + V^_ 

cxpj-y J dvdr' [qp+(v) - qp-(r) + a(r)]v(r - r') 
{qp + (r>)-qp_(r') + a(r>)} 
+i / dr [ip + (r)p + (r) + ^_(r)p_(r)] 



1 

iV+! 
1 

~nJ. 



^q 2 £ B v(0)/2 



dr 
dr 



J?(r)c 



Q(r)e 



3 -«/>+(r) 



-iijj- (r) 



N+ 



(3) 



The partition function becomes simpler upon transforma- 
tion to the grand-canonical ensemble according to 



E 



A 



N + +N- 




(4) 



where Ao is the bare fugacity which is related to the par- 
ticle chemical potential p, by Ao = e M . We assume here 
that the chemical potential of plus and minus ions is the 
same. Using the definition of the exponential function, 
e x = J2n=o xN /Nl, the grand-canonical partition func- 
tion can be written as 



with the field-theoretic action given by 

H[p+,p-> i>+^-] = 



(5) 



JL 
2 



3=1 



J drdr'[qp + (r) - qp-(r) + a(r))v(i 
[ q p + (r')-qp-(r')+*(r')] 
-i J dr [^+(r)/3+(r) +V-(r)p_(r)] 

-A / dr^(r) L-^+W +e -#-M 



(G) 
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where we defined a rescaled fugacity as A = X e q2tBV ^^ 2 / Af 
As follows from the definition of the grand-canonical par- 
tition function, Eq. (||) , the expectation value of the total 
particle number is given by 



(N+ + N-) = A 
A / dri?(r)(c 



dlnZ A 



A 



d\nZ x 
dX 



(7) 



As can be be seen, the expectation value of the particle 
number, (N), is independent of any multiplicative factor 
of the fugacity. The standard way of saying this is that 
the chemical potential is only defined up to a constant. 



2.1 Density description 

One type of approximations starts by eliminating the fluc- 
tuating fields tp+ and tp- from the problem. The simplest 
way of doing this is to replace the functional integral over 
these fields by the value of the integrand at the extremum 
or optimal value. The saddle point equations which deter- 
mine the optimal values of the fluctuating fields, i/j+ P an d 
ijj- P , are given by 



81>+(t) 







(8) 



and with a similar equation for ip_ . Using H as defined in 
Eq.(^|), the solutions are 



^(r)=thi(p + (r)/A)) 



(and a similar solution for ip_ ) which, upon insertion into 
the Hamiltonian Eq.(^j), leads to the saddle-point Hamil- 
tonian 

-f / drdr'[(7/7 + (r) - qp-(r) + a(r)]v(r - r') 
[ q p + (r')- q p_(r') + a(r')] 

dr p+(r) (ln[p+(r)/A]-l). 

+ / drp_(r)(ln[p_(r)/A]-l). (10) 



This expression is the starting point for local density func- 
tional theory Jll|. The extremum of Hamiltonian Eq. (|Io|) , 
determined by the equation 







(11) 



(and with a similar equation for pJ), is the mean-field or 
Poisson-Boltzmann density distribution. As is well-known, 
the mean-field equation ( |TTj ) can be solved in closed form 
only for a few cases, and for more complicated situations 
one has to rely on additional approximations. One pop- 
ular way of doing this is to minimize the Hamiltonian 



Eq.(p"o|) for a restricted functional form of the ion distri- 
bution, for example a box distribution. In the appendix 
Al we show how this can be done for a planar charged 
wall (where comparison with the closed-form mean-field 
solution is possible). The draw back of box models is that 
the ion distribution is discontinuous and thus rather ar- 
tificial and for more complicated charge distributions the 
box geometry becomes complicated. Another, more basic 
drawback of the density functional Eq.(10) is that corre- 
lation effects are not included on the local level. Non-local 
terms can in principle be taken into account by expanding 
around the saddle-point in Eq.(^|), but the resulting func- 
tional is complicated. In the next section we show that a 
variational approach is more useful when applied in the 
potential description. 



2.2 Potential description 

A more powerful variational formulation is possible by 
choosing a description in terms of the fluctuating poten- 
tials ip+ an d ip-. The starting point is the observation 
that the integral over the density fields p+ and p_ in Eq. 
(|5|) can be performed exactly. Introducing the new field 
4>{r) — [ijj + (r) — ip-(r)]/2, the partition function can be 
written exactly as 



v<t> 

Z,, 



(12) 



with 



(9) H[4>] 



dr 



[V<^(r)] 2 zcr(r)0(r) 



2XQ{v) cos0(r) 



(13) 

where we used the inverse of the Coulomb potential v 1 (r) = 
— V 2 £(r)/47r. The symbol Z v denotes the determinant of 
the Gaussian integral 

Zv = J 2?0cxp{-^- J drdr'^rKV-rW)} 

(14) 

and is a measure of the free energy of vacuum fluctuations. 

Next, we rescale all lengths by the Gouy-Chapman 
length p = 1/ (2Trq£B(J s ) according to r = pr , where a s de- 
notes the surface charge density of the charged object un- 
der consideration. We obtain the modified partition func- 
tion 

Z x = /^exp{-i/[0]/^} (15) 



with the rescaled action 



m - 1 1 



i(V0(f)) 2 +2<? 



(f)a(f) 



AD{r) cos 0(f) 



where we used the rescaled charge distribution a(r/ p) = 
pa(r)/a s and ion exclusion function Q[r/ p) = Q(r). The 
rescaled fugacity A is defined by 



A = 8nXp 3 S 



2A 



(17) 
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and is determined by the normalization condition (in the 
case of fixed particle number) 



(JV 4 
A 

a-kE 



N-) 



A 



dlnZ A 
dA 



— / df 17(f) (cos 0(f)). 



(18) 



As one can see, there are two parameters left in the field- 
theoretic formulation, namely the coupling parameter 



E = 27re 2 B cr s q 3 



(19) 



which measures to which extent deviations from the saddle- 
point are important, and the ion fugacity A. Infinitely far 
from any charged object the potential is constant and the 
concentration of minus and plus ions is equal and given 

* N - A ■ *) (20, 



V 



V 



which, using the definition of the screening length, 



(8nq 2 £BC s ) x l 2 and the Gouy Chapman length can be 
rewritten as 

■ 2 = k 2 = A(cos(f>). (21) 



K 2 /J,~ 



The rescaled fugacity A thus measures the ratio between 
the screening length and the Gouy-Chapman length: For 
k 3> 1 the screening by salt is dominant, non-linear ef- 
fects are unimportant and Debye-Huckcl theory is valid, 
for k< 1 on the other hand non-linear effects are impor- 
tant and Debye-Hiickel theory fails. It is the latter case 
where our variational theory becomes important. 



3 Variational Theory 

The action Eq.(|l|) is too complicated to be solved exactly. 
Previously, perturbative theories have been proposed for 
the weak-coupling (small E) regime in terms of a loop- 
expansion around the saddle-point solution (which corre- 
sponds to the solution of the non-linear Poisson-Boltzmann 
equation) jl4j and for the strong-coupling (large E) regime 
in terms of a virial expansion |lE|, |l6|, [l7t . Both expansions 
are of limited usefulness because of their bad convergence 
properties as one moves into the intermediate regime be- 
tween the two asymptotic limits j^]. It is this point where 
variational theories come in as a robust means to obtain 
reliable, useful results. 

The standard Gibbs variational procedure consists in 
minimizing the following variational free energy 



Fdbbs — J~o — (Ho — H)o/^ 



(22) 



with respect to the variational parameters included in the 
variational Hamiltonian Hq. As discussed in detail in Ap- 
pendix B the Gibbs variational principle is equivalent to 
the first-order perturbational variation on the observable 
which is conjugate to the variational parameter. The per- 
turbational variation scheme can be systematically im- 
proved by going to higher order, the Gibbs variation can- 
not. Since we stay at the first-order level throughout this 



paper, we will use the Gibbs variational notation, though 
one should keep in mind that this can be interpreted as 
the first order result of a perturbational variational calcu- 
lation. 



3.1 Most general Gaussian variation 

The most general, yet exactly tractable variational Hamil- 
tonian, is a Gaussian one 

H [4>\ = \ J [0(r)- i 0o(r)H- 1 (r,rO[0(rO-«0o(r')] 

(23) 

where the Gaussian kernel v and the mean potential 
tpo are variational parameters. The variational free energy 
according to Eq. ( f22| ) reads 

SGibbs = --Trlogvo - / 8?r „ 



S{r - r 

r' 

<?(r)<M r ) 



, V r Vr'Uo(r,r') 



8nE 



2nE 



A 
AttE 



q e -s« (p,r)/a co8h ^j < ( 24 ) 



The variational equations are given by SJ-aibbs/S^oir) = 
and 5J r Gibbs/5vQ 1 {r, r') = and read 



= + —n e -=Mr,r)/2 sinh0o (25) 

4-7T 2-7T An 



V 2 z; (r,r') A _ s 



5(r-r') 



+ _ e -3vo{r,v)/2 cogh ^ Q ( j/) 

Att 



(26) 



These self-consistent equations are quite complicated 
even for the simplest charge distribution. In the following 
we will devise a simpler, one-parameter variational scheme 
which will lead -as a byproduct- to a very simple and 
clear definition of the renormalizcd or effective charge of 
an object. 



3.2 Variational Charge Renormalization 

As we will show in the following, non-linear effects are 
captured in an essentially exact way already by a much 
simpler variational form, which as the only variational 
parameter contains a charge-renormalization factor. The 
variational Hamiltonian reads 



Z7T 



i(V0(f)) 2 +t?? (f)0(f)a(f) + ^M 



( 2? ) 

where the variational function rj (r) multiplies the charge 
distribution and therefore corresponds to a spatially de- 
pendent charge renormalization factor. The Hamiltonian 
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is quadratic in the fluctuating field and therefore all expec- 3.3 Homogeneous variational approach 

tation values can be calculated by using Green's function 

formalism. Also, the geometry-dependent factor J?(r) is The homogeneous variational approach is obtained by choos- 

for the present formulation omitted and the salt concen- ing a charge renormalization which is uniform in space, i.e. 

tration dependent term proportional to k 2 has the same rj(r) = rj . We redefine the bare potential as 
strength throughout the space. This is certainly a simpli- 
fication, which however allows to use a simple translation- , , . 1 f _ , , , , 
ally invariant Green's function (future elaborations should 2n J, *" ) Vjm y r T ) \ > 
probe variational forms with a spatially dependent mass 

term). Using the Debye-Hiickel potential and the Hamiltonian can be rewritten as 

_ If 

v DH (r) = '—- (28) ^o[(/>] = ~ [Hr) + ir]oM r )] v Bii( F ^')[^( r ') + l 'noM F ')} 

and the potential distribution +— / a(r)<j) (r). (35) 

47T J r 

M r ) = ^ J ?7(r')o-( r ') u DH(r - r') (29) The variational free energy according to Eq.(||) reads 

the variational Hamiltonian Eq.@ can be rewritten as F&bbs = ^dh + — ^= 1(2- r] Q )r] a(r)(po(r) (36) 

4?r^ J t 

M4>] =\ [0(r) + i o (r)]t- D i (r,r')[0(r') +^o(r')] [ fog.tfftr) - S«dh(0)1 

+ JZ I v(*)Z(?)<h( r )- ( 3 °) --^ / ^(r)cosh[ry o 0o(r)] 



The variational free energy according to Eq.fl22J) is easily The var i a ti nal equation follows by taking a derivative 
calculated and reads with respect to the charge renormalization factor, dJ^cabs I dvo 

0, and is given by 



FGibbs -^dh + ^ j(2 - f/(r))?(r)^o(r) (31) 



~ K 2 , 2(1-%) / cr(r)to(r) (37) 



K 



8irE; 



S(r)-5WDH(0)] 



- / f?(r)cosh[0 o (r)] 



« 2 / 0o (r) ^(r)sinh[r/ o 0o(r)] - r] ^ (i 



In the following we will solve this self-consistent equa- 
where we used the short-hand notation tion for the charge-renormalization factor 770 in the planar, 



•^DH = - hi 



T>6 e~f*S *(»>dh('-0*('')/ 2S 



cylindrical, and spherical geometries. 



1 f d 3 q ,, 99n , N 3.3.1 Planar case 

' Mog(q 2 + K 2 ) (32) 



2 J (2tt) 3 

For the calculation of the uniformly charged planar wall 
and fixed the fugacity A according to Eq.((2|). The vari- we assume that salt is present on both sidesof the wall, 
ational equation follows by taking a functional derivative the geometry factor therefore equals unity Q(r) = 1 (it 
of the free energy with respect to the charge renormaliza- is an easy exercise to generalize the calculation to the 
tion function, SFGibbs / ^vi*) — 0; ancl reads after a few case where salt is present only on one side of the charged 
algebraic manipulations wall). With the bare charge distribution <r(r) = a s 5(z) the 

rescaled charge distribution becomes cr(r) = S(z). The 

fl - n(r))a(r)6n(r) = bare P otential is according toEq.(§|) given by (f) (z) = 

v /v v JV0{ ' e~ KZ /k. The integrals in Eq.@ can all be performed ex- 

z.2 p actly, and the resulting equation becomes 
— J o (r) (!7(r)sinh[0 o (r)] - (fo(r)J . (33) 

1 - ^ - — (cosh[%/£] - 1) = 0. (38) 
This integral equation is to be solved by a suitably cho- 2 770 

sen function ?7(r). Since in this paper we will be only con- , . 

cerned with situation where translational invariance holds, In the llmlt K > 1 the solutlon 1S 
we will in the following section specialize the variational \ 
approach to the homogeneous case. r /o — 1 — 24^2 
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and therefore only weak charge-regulation takes place, as 
expected for a weakly charged surface or high salt concen- 
trations. In the opposite limit of a highly charged surface 
or at low salt concentrations, for k <C 1, the solution is 



-k ln(ft) 



(40) 



These results can be directly compared with the Poisson- 
Boltzmann results derived explicitly in Appendix A2. The 
quantity one wants to compare is the surface potential, 
since this is conjugate to the surface charge density, i.e. the 
variational parameter, and one can reinterpret the Gibbs 
variational principle as a first-order perturbational varia- 
tion on the expectation value of the surface potential (as 
is explained in more detail in Appendix B). Within our 
scheme, the surface potential is defined by 

V^o = mM* = 0) 

and is given in units of fc^T and valency q. For the pla- 
nar case, we obtain -0o = Vo/k, and therefore find for the 
surface potential r/o/k = 1/k at high salt (k 1) and 
Vo/k — ~ at low salt (R -C 1). These are the surface 
potentials one obtains from the full non-linear Poisson- 
Boltzmann theory (except a factor of two for the low- 
salt case) . This shows that our variational procedure cap- 
tures non-linear effects accurately and justifies to con- 
sider more complicated geometries, where the Poisson- 
Boltzmann equation can only be solved numerically. 



3.3.2 Cylindrical case 

We assume a charged cylinder with radius R and linear 
charge density r. The surface charge density follows as 
a s — t/(2ttR). Defining the Gouy-Chapman length (like 
in the planar case) as fx = l/(2i:qa s ), we obtain /i = 
R/(qr£B)- The rescaled cylinder radius, R = R/ ' [i = qr£s 
is thus equal to the Manning parameter ||: For R > 1 one 
expects counterions to be bound to the cylinder even in 
the limit of vanishing salt concentration, for R < 1 no 
counterions stay bound in that limit. In the following, we 
assume the cylinder to be penetrable for ions, the geom- 
etry function is f2(r ) = 1 for all values of f. Using cylin- 
drical coordinates, the potential 0o as defined in Eq.(p4[) 
satisfies the differential equation 

1 c) c) 

--jpr—Mr) + 25(f -R) + k 2 M?) = 
r Or Or 

with the boundary conditions (f>o(R — 0) = <f>o(R + 0) and 
(j>' (R + 0) — (f>' (R — 0) = 2. The complete solution is con- 
structed from the two solutions of the homogeneous dif- 
ferential equation, <po{r) — Ko(kr) and <po(f) — Io(kf), 
and reads 

M?) = 2RI Q (kR)K (kf) (41) 



for f > R and 



for f < R. The solution of the integral equation ( p7| ) can be 
easily obtained graphically by calculating r]Q as a function 
of kR and rjo/k via numerical integration followed by a 
parametric plot. In Fig. la we show the renormalization 
factor r)o as a function of the Manning parameter R for 
fixed values of the rescaled inverse screening length kR = 
5,1,0.2,0.01 (thin solid lines from top to bottom). The 
thick solid line denotes the Manning limit 



77o = 1/R = l/(r£ B q) 



(43) 



which is very slowly approached as the salt concentration 
decreases 0. Indeed, an asymptotic analysis of the self- 
consistent equation (37) with the cylinder potential Eqs. 
©II) gives for kR < 1 



1 I \nR 

^il 1+ 21n(l/ K i?) 



(44) 



which shows that the approach to the Manning limit is 
logarithmically slow. In the unrescalcd units, the effective 
cylinder charge density is 



T e ff = TT]q 



1 



1 



In risq 
2\n(l/ kR) 



(45) 



In Fig. lb we show the data for 770 for fixed Manning pa- 
rameter R = 3 (which roughly corresponds to a fully 
charged vinyl-based polyelectrolyte with a distance of 0.254 
nm between charges along the backbone) as a function of 
kR together with the asymptotic prediction Eq.(|4^) as a 
broken line in a semi- logarithmic plot and in Fig.lc as 
a function of — l/ln(«;i?). It is seen that the asymptotic 
law is valid for low salt concentrations satisfying kR < 0.1. 
For large salt concentrations, asymptotic analysis predicts 
that the charge-renormalization vanishes according to 



1 



1 



24k 2 



(46) 



Mr) = 2RK (kR)I {kr) 



(42) 



(identical to the result for the planar case) which is shown 
in Fig. lb as a dotted line. 



3.3.3 Nonlinear electrostatic persistence length 

The electrostatic contribution to the bending rigidity has 
been calculated on the linear level |],|l(| but also on the 
non-linear level using the full Poisson-Boltzmann equa- 
tion for a bent cylinder[jl9[|2^]. Here we re- address the 
relevance of non-linear effects for the persistence length 
of highly charged polyelectrolytes. In Fig. Id we plot the 
charge renormalization factor on double logarithmic scales. 
Over roughly a decade in kR, which corresponds to two 
decades in the salt concentration, the charge renormaliza- 
tion factor can be approximately described by a power law 
770 ~ n a with a « 0.3 (shown as a solid line). This range 
of salt concentrations is quite relevant for experiments Jl2[ , 
since it has been reported that the experimentally mea- 
sured salt-concentration dependence of the persistence length 
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-l/log(KR) 



1 10 

kR 



Fig. 1. a) Charge renormalization factor rjo for a charged cylinder of radius J? as a function of the Manning parameter R = qtBT 
for fixed values of the rescaled inverse screening length kR — 5, 1, 0.2, 0.01 (thin solid lines from top to bottom). The thick solid 
line denotes the Manning limit Eq.(|43|). b-d) Results for fixed Manning parameter R — 3 as a function of the rescaled radius 
compared with the asymptotic low-salt prediction Eq.(|44[) (broken line) and the asymptotic high-salt prediction Eq.(^) (dotted 
line). The solid line in d) corresponds to the power-law 770 ~ k o z . 



deviates from the Odijk-Skolnick-Fixman prediction, ac- 
cording to which 

iosF = -fp- (47) 

with P — 2|],[l(|. Replacing the bare line charge density 
r in Eq.fi^) by the effective charge density 



T e ff = ?7ot 



(48) 



(which is permissible since to leading order in an expan- 
sion for large bending radii modifications of the charge 
renormalization due to polymer bending are irrelevant) 
one obtains a modified exponent /3 = 2 — 2a 1.4 which 
is closer to the experimentally measured exponents]!^. 
Other effects that have been shown to reduce the expo- 
nent (3 are relevant for weakly charge chains which are 
crumpled at small length scales) 11 2Up2],23 although re- 



cent simulations using linear Debye-Huckel potentials [24, 
|25| as well as general Gaussian variational calculations! 26 
tend to confirm the original Odijk-Skolnick-Fixman pre- 
diction with P = 2. It might therefore turn out that the 



true explanation for the discrepancy between experimen- 
tal and theoretical predictions for the electrostatic persis- 
tence length is due to non-linear effects which had been 
neglected in most theoretical treatments but are of course 
omnipresent in experiments. Needless to say, we expect 
modification of a whole number of other scaling relations 
for polyelectrolytes which depend on the line charge den- 
sity r. 



3.3.4 Spherical Case 

We now treat a sphere of charge Z and radius R, leading 
to a surface charge density a s — Z/(AnR 2 ) and Gouy- 
Chapman length /i = l/(2Trq£B<7 s )- The rescaled radius is 
given by R = R/y, = ql B Z/{2R). 

Using spherical coordinates, the potential <po as defined 
in Eq. (|34j) satisfies the differential equation 



f 2 dr df 



2 d -M?) + 2<*(r -R) + K 2 Mr) = 
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with the boundary conditions <j>o(R — 0) = 4>o(R + 0) and 
<j)' {R + 0) - (j)' {R - 0) = 2. The two solutions of the ho- 
mogeneous differential equation are 4>o{f) — e Kr /f and 
cf>o (f) = c~ Rr /r, from which we obtain the solutions 

Mr) = 2flsinh(/t.R)e- sf /(fir) (49) 

for r > R and 

0o(f) = 2Re~ iiR smh(k?)/(k?) (50) 

for f < R where we have again assumed that salt ions are 
present everywhere, i.e. f2(r) = 1 throughout space. 

In Fig. 2a we show parametric plots of the charge renor- 
malization factor for kR = 100, 10, 1, 0.4 (from top to bot- 
tom) and in Fig.2b for kR = 0.4, , 0.1, 0.01, 0.001 (from 
bottom to top). In contrast to the cylindrical case, charge 
renormalization is most pronounced at intermediate salt 
concentrations corresponding to a screening length which 
roughly equals the sphere radius For large salt con- 
centrations, asymptotic analysis predicts that the charge 
renormalization vanishes according to 



in agreement with the planar and cylindrical geometries. 
For low salt concentrations, an asymptotic analysis of the 
self-consistent equation ([37]) with the spherical potential 
Eqs. ©HI gives for kR < 1 

77 ~-^lnfJ— ). (52) 

2R \Ray y ' 

In unrescaled units this leads to an effective sphere charge 

Zcs = v ° z ^^L in {w§) (53) 

which only depends very weakly on the bare sphere charge 
(in agreement with previous qualitative arguments^). In 
Fig. 2c we compare data for a sphere with rescaled radius 
R = Rj \i = 10 with the asymptotic low-salt prediction 
Eq.(p2|) (broken line) and the high-salt prediction Eq.(|5l|) 
(dotted line). Note that the asymptotic low-salt predic- 
tion is only valid in a finite salt-concentration limit: for 
very low salt concentration, the counterions are entrop- 
ically driven away from the sphere and the bare charge 
is recovered. This might be detectable in light-scattering 
experiments on de-salted charged-colloid suspensions. 



4 Conclusions 




Fig. 2. Results for the charge renormalization factor rjo 
for a sphere for rescaled inverse screening lengths a) kR — 
100,10,1,0.4 (from top to bottom) and b) for kR = 
0.4, ,0.1, 0.01, 0.001 (from bottom to top), c) Charge renor- 
malization factor for fixed R = R/fi = 10 together with the 
asymptotic prediction Eq.(^) (broken line) and the prediction 
Eq.@ (dotted line). 



We have introduced the framework for a general varia- 
tional technique which allows to capture non-linear ef- 
fects of counterion distributions around highly charged 
objects in an essentially exact way. The scheme builds 
on the field-theoretic formulation of fluctuating counteri- 
ons and uses the effective surface charge as a variational 
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parameter. Since we use a Gaussian variational Hamilto- 
nian, the resulting interactions between charged particles 
are of the Debye-Hiickel type, but with the bare charges 
replaced by effective ones. This scheme is thus very close 
to the treatment by Kjellander et al.^Tj, where a similar 
procedure was used within an integral-equation scheme. 
We show how in principle to systematically improve the 
leading term results by using perturbational variational 
methods. We tested our approach for the planar geometry, 
where comparison with the analytical solution of the non- 
linear Poisson-Boltzmann equation is possible, and then 
presented a detailed analysis of the cylindrical and spheri- 
cal case. In the future, we will apply these methods to the 
interaction between charged bodies (including the experi- 
mentally relevant case of a charged cylinder adsorbing on a 
charge plane Q), soft (fluctuating) charged polymers and 
membranes. We also plan to investigate higher-order per- 
turbation variational terms which will allow to include ef- 
fects at high coupling parameter S where deviations from 
the mean- field like behavior are important. 



A Appendix: Charged planar wall 

In this section we will treat the single charged wall im- 
mersed in a salt solution using a box model and using the 
Poisson Boltzmann equation. 



A.l Box model 



The excess entropy of the layer of condensed counterions 
is 



(J.s 



S <l < \ p° + + —\ \ lot; 



(57) 

Minimization of the free energy F = U — S leads for high 
salt concentrations (for up — k 1) to the optimal layer 
thickness 



An£ B q 2 p 



2i 



6k 



(58) 



which is proportional to the screening length k 1 (as ex- 
pected) and for low salt concentrations (for up <C 1) to 



2ir£ B qa< 



= 3p 



(59) 



and thus to a layer thickness which differs from the Gouy 
Chapman length p only by a numerical factor. The surface 
free energies are given by 



iqnp 



for high salt and 



for low salt. 



q \np 



(60) 



(61) 



The mean-field Hamiltonian in terms of the ionic densities, 
Eq.(|l0|), can be used as a starting point for a simplified 
variational treatment using a box model, where the true 
ionic density distribution is replaced by a step profile. The 
fugacity can be fixed by the ionic density in the bulk where 
we assume the electrostatic interactions to vanish. Mini- 
mization Eq.([Io|) with respect to the density of plus ions 
yields 

A FT 

— = H P+ /X) = (54) 
Sp + 

and thus p + = A. We conclude that the fugacity is given 
by the ionic density in the bulk reservoir, i.e. A = p h + . 

As a box model we assume a negatively charged planar 
wall of surface charge density cr s which is neutralized by a 
uniformly charged box of width d. We assume that the wall 
charge is neutralized by an excess of positive ions, i.e. the 
density of negative ions stays constant. The electric field 
distribution thus is 



eE d - z 

= An£ B a s — j— 



(55) 



and the electrostatic energy (per unit area and per k B T) 
follows as 



U 



dzE 2 



1 



2k B T J 8n£ B 
= 2n£ B a 2 d/3. 



dz 



eE 



(56) 



A. 2 Poisson-Boltzmann equation 

The Poisson-Boltzmann equation follows from minimizing 
the Hamiltonian Eq. ( |l6| ) with respect to the potential dis- 
tribution. Defining the real electrostatic potential ip = t<f> 
the saddle point equation can be written as 



ip"(z) + 2S(z) - k 2 sinh V> = 0. 
The solution is 



rj>(z) = 2 log 



1 — 7e 



1 



■ 7e 



(62) 



(63) 



where the constant 7 is determined by the boundary con- 
ditions. In the case when the charged wall is impenetrable 
to ions, the boundary condition is ip'(0) = —2, which leads 
to the equation 7 2 — 1 = 27ft. For low salt concentrations 
the constant follows as 7 ~ — 1 + ft and the surface poten- 
tial is given by 

Vo-21n(l/ft). (64) 

The surface free energy in this limit is F = tpoa s /(2q) = 
(a s /q) ln(l/ft) and thus differs from the box-model calcu- 
lation only by a prefactor. For high salt concentrations the 
constant is 7 = — I /(2k) and the surface potential is 



ip — 2/ft. 



(65) 
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free energy, T\ > T . The approximate free energy is there- 
fore given by minimizing T\ with respect to the set of 
variational parameters, 

-^Gibbs = minJ r i(x). (73) 

X 

However, it is not easy to systematically improve this re- 
sult (since the property of the perturbative free energy 
to be bound by the true free energy is lost at higher or- 
der), other than by choosing a more refined variational 
Hamiltonian, which is often not possible. The perturba- 
tional variation, which has been introduced by Edwards 
for the calculation of the scaling behavior of self-avoiding 
polymers p9|, consists in requiring the perturbative free 
energy of any order to be equal to the zeroth order free 
energy, 

.Fo(x) = jr„(x) (74) 
which is equivalent to the condition 

n . m 

E ~T ( Ym )o = 0- (75) 
* — ' m! 

m— 1 

This equation has to be solved by fixing the set of varia- 
tional parameters x. To do this, it is convenient to expand 
the parameters in a power series x = Y17=i C~ lx ^ an d to 
solve Eq. (|75|) term by term. (If the variational parameter 
is not expanded in a series but Eq.(f75|) is solved directly, 
one recovers the naive perturbation result. The hope is 
therefore that the expansion of the variational parameters 
improves the convergence behavior of the series for the 
free energy.) One notes that the prediction for the free 
energy from the Gibbs variational approach, Eq.(f73|), and 
the perturbational variation, Eq.([73), are not equal, not 
even at the leading order. It is a priori not clear which 
approximation closer to the exact free energy. 

Often one is not interested in the free energy itself but 
in a certain observable. We define an arbitrary observable 
Xj as being conjugate to the field Xj according to 

The first few terms of the expansion of the expectation 
value of Xj in powers of £ are 
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The surface free energy is F = a s /(kq) and is also repro- 
duced by the box-model calculation. In the case when the 
surface is penetrable to ions, as is used in our variational 
calculations, the surface potential in the low salt limit is 
not modified to leading order, 

^o~21n(l/«) 1 (66) 

while in the high salt limit one obtains 

</>o 1/k. (67) 

B Appendix: Connection between the Gibbs 
Variational Principle and Perturbative 
Variation 

Let us start by considering the free energy 

T = - In J V<j)e- n (68) 

where 7i depends on the fluctuating field cj). We intro- 
duce a variational Hamiltonian 7io( x ) which depends on 
the fluctuating field 4> and on the vector x, the compo- 
nents of which may act as variational parameters and/or 
as generating fields to calculate observables. Adding and 
subtracting the variational Hamiltonian to the action, the 
free energy can be written as 

j (69) 

where we defined 

Y = Hq{x)-H (70) 

and introduced a dummy parameter £ which is used to 
count orders in a perturbation series (and will eventually 
be set to unity). The free energy can now be expanded in 
cumulants as 

71 /~rn 

T n = T a {*)-Y,—(Y m )» (71) 
'-^ ml 

m—l 

where .Fo(x) denotes the free energy of the variational 
Hamiltonian 

^o(x) = - In J P0e^ Ho(x) (72) 

and (. . .) denotes the cumulant average with respect to 
the variational Hamiltonian. Clearly, the exact free en- 
ergy is approached in the infinite limit of the series, i.e. 
T = Tn^ao. The standard perturbation technique is of- 
ten slowly converging for problems of interest. Variational 
techniques typically perform much better. The standard 
Gibbs- variational technique exploits the fact that the first 
order perturbation result T\ is always larger than the true 



(Xj) = <X,) + C [<X,y> - (Xj)o(Y)o] (77) 
+ y [(XjY^o - 2{X j Y) {Y) - (XMY 2 )o + 2{X j ) {Y)*] . 



As before, the perturbative variation consists of demand- 
ing that the corrections to the leading contribution of the 
expectation value vanish, i.e. 

(Xj) = (Xj) , (78) 

which can be solved by adjusting one or more of the vari- 
ational parameters x. 

A quite instructive connection exists between the Gibbs 
variational method and the perturbative variation: In the 
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case when the Gibbs variational parameter is the field 
Xj conjugate to the observable Xj, the Gibbs variational 
equation and the first order perturbative variation Eq.(f78|) 
are identical. This means that the higher order perturba- 
tive variational method on the expectation value can be 
viewed as a generalization of the Gibbs variational tech- 
nique. One might hope to approximate the true free energy 
better, though the true free energy ceases to be a strict 
lower bound for the perturbational free energy expression. 
Some comments are in order: 

i) In principle, the calculated observable in equation 
) does not need to be conjugate to the variational pa- 
rameter. In any case, it is clear that in principle the observ- 
able can be calculated to any wanted precision by taking 
enough terms in the perturbation- variational treatment. 

ii) Several expectation values can be calculated, for 
which an equal number of variational parameters are needed 
to solve the uniform expansion condition. 

v) It is not clear that equation ([78|) always has a solu- 
tion. In that case one could try to find the best estimate 
to the solution. 
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